Demographics and reproductive biology of H. schistosus from bycatch on the west coast of India
Analysis for Manuscript
Demographics of H. schistosus across time
Sample size
The number of snakes (N) sampled from bycatch of trawlers, shore seines and gillnetss.
snakes%>%
group_by(Species)%>%
count(name = "N")| Species | N |
|---|---|
| Hydrophis curtus | 282 |
| Hydrophis schistosus | 967 |
The number of trips (n) and fishing effort (haul-hours) sampled for sea snakes in the current study.
snakes%>%
filter(Gear.Type != "")%>%
select(Date, Gear.Type, No..of.Hauls, Average.Haul.Duration..Hours., Tow.duration.hours)%>%
distinct()%>%
# Calculating tow duration
mutate(Tow.duration.hours = case_when(!is.na(Tow.duration.hours) ~ Tow.duration.hours,
is.na(Tow.duration.hours) ~ No..of.Hauls*Average.Haul.Duration..Hours.))%>%
group_by(Gear.Type)%>%
summarise(n = n(), # counting number of trips sampled
haul.hours = sum(Tow.duration.hours, na.rm = T)) # total fishing effort| Gear.Type | n | haul.hours |
|---|---|---|
| GillNet | 340 | 535.9633 |
| Rampan | 46 | 190.4500 |
| Trawler | 76 | 285.1600 |
Age structure
# calculating the mean SVL across years
age.yr <- snakes%>%
group_by(Species, Year)%>%
summarise(N = n(),
mean = mean(Snout.to.Vent..cm., na.rm = T))
# SVL cut off for age classes
maturtity <- snakes%>%
group_by(Species)%>%
count()%>%
mutate(juv = 35,
adult = ifelse(Species == "Hydrophis curtus", 54, 65))Distribution fo SVL across years
# plotting distribution of SVL across years
snakes%>%
filter(Snout.to.Vent..cm. > 20)%>%
ggplot(aes(Year, Snout.to.Vent..cm.))+
geom_violin(fill = "grey")+
geom_boxplot(width = 0.1)+
geom_hline(data = maturtity, aes(yintercept = adult), linetype = "dashed")+
geom_hline(data = maturtity, aes(yintercept = juv), linetype = "dotted")+
#stat_summary(fun.data = "mean_sdl", geom = "pointrange", size = 1)+
geom_label(data = age.yr, aes(Year, 10, label = N))+
facet_wrap(~Species, ncol = 1, scales = "free_y")+
labs(y = "Snout to vent length (cm)")+
theme(strip.text = element_text(face = "italic"))ggsave(last_plot(), filename = "./Figures/figure1.tiff", height = 6, width = 8)We did not find any H. curtus neonates in our sampling. Majority of H. curtus were juveniles. H. schistosus populations consisted of mostly adults.
Proportion of developmental classes across years
pop.yr <- snakes%>%
filter(!is.na(Class))%>% # fix 2016 svl
group_by(Species, Year)%>%
# Counting number of individuals in each developmental class
count(Class)%>%
complete(Class, fill = list(n = 0))%>%
# Total number sampled in each year
mutate(N = sum(n))%>%
group_by(Year, Species, Class)%>%
# Proportion of each developmental class
mutate(prop.age = n/N)
# Summary stats
pop.yr%>%
group_by(Species, Class)%>%
skimr::skim(prop.age)%>%
skimr::yank("numeric")%>%
select(c(Species, Class, mean, sd))Variable type: numeric
| Species | Class | mean | sd |
|---|---|---|---|
| Hydrophis curtus | Adult | 0.43 | 0.31 |
| Hydrophis curtus | Juvenile | 0.57 | 0.31 |
| Hydrophis curtus | Neonate | 0.02 | NA |
| Hydrophis schistosus | Adult | 0.88 | 0.09 |
| Hydrophis schistosus | Juvenile | 0.14 | 0.07 |
| Hydrophis schistosus | Neonate | 0.03 | 0.02 |
Testing change in proportion of developmental classes across years
pop.yr%>%
filter(Year != 2019)%>%
group_by(Species, Class)%>%
nest()%>%
mutate(test = map(data, ~prop.test(.$n, .$N)), # Z - test for proportion
sumr = map(test, broom::tidy))%>%
select(sumr)%>%
unnest()%>%
select(Class:p.value)| Species | Class | estimate1 | estimate2 | estimate3 | statistic | p.value |
|---|---|---|---|---|---|---|
| Hydrophis curtus | Adult | 0.8750000 | 0.3125000 | 0.3551402 | 9.2546983 | 0.0097807 |
| Hydrophis curtus | Juvenile | 0.1250000 | 0.6875000 | 0.6261682 | 8.8509715 | 0.0119684 |
| Hydrophis curtus | Neonate | NA | NA | NA | 97.2336449 | 0.0000000 |
| Hydrophis schistosus | Adult | 0.8846154 | 0.7894737 | 0.8613445 | 5.9949721 | 0.0499124 |
| Hydrophis schistosus | Juvenile | 0.0769231 | 0.2105263 | 0.1239496 | 9.7219255 | 0.0077430 |
| Hydrophis schistosus | Neonate | 0.0384615 | 0.0147059 | NA | 0.0189777 | 0.8904305 |
Age structure does not change significantly over a four year period from 2016 to 2019.
Change in SVL distribution over years
library(car)
snakes%>%
filter(Year != "2019")%>%
group_by(Species)%>%
select(Year, Snout.to.Vent..cm.)%>%
nest()%>%
mutate(mod = map(data, ~lm(Snout.to.Vent..cm. ~ Year, data = .)),
sumr = map(mod, broom::tidy),
stat = map(mod, car::Anova))%>%
select(stat)%>%
unnest()| Species | Sum Sq | Df | F value | Pr(>F) |
|---|---|---|---|---|
| Hydrophis schistosus | 9866.294 | 2 | 21.86203 | 0.00e+00 |
| Hydrophis schistosus | 168785.548 | 748 | NA | NA |
| Hydrophis curtus | 2641.334 | 2 | 10.07170 | 7.05e-05 |
| Hydrophis curtus | 24258.408 | 185 | NA | NA |
There was a significant change in SVL distribution of species across years.
month.svl <- HS%>%
filter(Year == "2018", # data not sufficient for other years
Snout.to.Vent..cm. > 20)%>% # removing erroneous data
mutate(Month = factor(Month, levels = month.name))%>%
complete(Month)%>%
group_by(Month)%>%
# Calculating mean SVL in each month
summarise(mean.SVL = mean(Snout.to.Vent..cm., na.rm = T))
# Month of observed births
births <- data.frame(Species = "Hydrophis schistosus", Month = "April")# plotting distribution of SVL across months
HS%>%
filter(Year == "2018", # data not sufficient for other years
Snout.to.Vent..cm. > 20)%>% # removing erroneous data
mutate(Month = factor(Month, levels = month.name))%>%
complete(nesting(Species), Month)%>%
droplevels()%>%
ggplot(aes(Month, Snout.to.Vent..cm.))+
geom_violin(fill = "light grey")+
#geom_point(data = month.svl, aes(x = Month, y = mean.SVL), size = 3)+
geom_boxplot(width = 0.1)+
geom_segment(data = births, aes(x = Month, xend = Month, y = 0 , yend = 10), #marking births
arrow = arrow(length = unit(0.25, "cm"), ends = "first"), size = 1)+
geom_text(data = births, aes(x = Month, y = 20, label = paste("Observed \n births")))+
geom_vline(aes(xintercept = "June"), size = 1)+#start of the monsoon ban
geom_vline(aes(xintercept = "August"), size = 1)+#end of the monsoon ban
geom_text(aes(x = "July", y = 80, label = "Monsoon ban"))+
scale_x_discrete(guide = guide_axis(n.dodge = 2))+
labs(y = "Snout to vent length (cm)")Proportion of neonates increases in March to May. Birth and neonates were observed around the same time.
Proportion of jeuveniles in each month of 2018
snakes%>%
filter(Year == "2018", # data not sufficient for other years
Snout.to.Vent..cm. > 20)%>% # removing erroneous data
mutate(Month = factor(Month, levels = month.name))%>%
group_by(Species, Month)%>%
summarise(prop.neonate = sum(Snout.to.Vent..cm. < 40)/n())%>%
spread(Month, prop.neonate)| Species | January | February | March | April | May | October | November | December |
|---|---|---|---|---|---|---|---|---|
| Hydrophis curtus | 0 | 0 | 0.50 | NA | NA | 0 | 0.0769231 | 0 |
| Hydrophis schistosus | 0 | 0 | 0.02 | 0.05 | 0.0555556 | 0 | 0.0000000 | 0 |
Sex ratios
Plotting sex ratios for each species across years
snakes%>%
filter(Sex == "Male" | Sex == "Female")%>%
group_by( Species, Year)%>%
summarise(N = n(),
females = sum(Sex == "Female"),
prop.female = females/N)%>%
ggplot(aes(Year, prop.female, fill = Species))+
geom_col(width = 0.5, col = "black", position = position_dodge())Proportion of females to entire population stays constant over the sampling period.
Proportion of females encountered by species
snakes%>%
group_by(Species)%>%
filter(Sex == "Male" | Sex == "Female")%>%
summarise(females = sum(Sex == "Female"),
N = n())| Species | females | N |
|---|---|---|
| Hydrophis curtus | 102 | 212 |
| Hydrophis schistosus | 304 | 618 |
Testing shift in sex ratios of sea snakes in bycatch
# is proportion of females different from 0.5?
snakes%>%
filter(Sex == "Male" | Sex == "Female")%>%
group_by(Species)%>%
summarise(females = sum(Sex == "Female"),
N = n())%>%
group_by(Species)%>%
nest()%>%
mutate(test = map(data, ~prop.test(.$females, .$N, p = 0.5)), # Z - test for proportion
sumr = map(test, broom::tidy))%>%
select(sumr)%>%
unnest()%>%
select(-c(method, alternative))| Species | estimate | statistic | p.value | parameter | conf.low | conf.high |
|---|---|---|---|---|---|---|
| Hydrophis curtus | 0.4811321 | 0.2311321 | 0.6306857 | 1 | 0.4125066 | 0.5504525 |
| Hydrophis schistosus | 0.4919094 | 0.1310680 | 0.7173273 | 1 | 0.4518628 | 0.5320580 |
Sex ratio is not different from 0.5; p = 0.71.
Mortality in bycatch
Overall mortality rate
snakes%>%
group_by(Species)%>%
count(Condition.at.encounter..D.A.)%>%
# Calculating mortality rate
mutate(N = sum(n),
prop.dead = n/N)%>%
filter(Condition.at.encounter..D.A. == "Dead")%>%
select(-Condition.at.encounter..D.A.)| Species | n | N | prop.dead |
|---|---|---|---|
| Hydrophis curtus | 120 | 282 | 0.4255319 |
| Hydrophis schistosus | 174 | 967 | 0.1799380 |
Testing difference in mortality across species
snakes%>%
group_by(Species)%>%
count(Condition.at.encounter..D.A.)%>%
# Calculating mortality rate
mutate(N = sum(n),
prop.dead = n/N)%>%
filter(Condition.at.encounter..D.A. == "Dead")%>%
select(-Condition.at.encounter..D.A.)%>%
ungroup()%>%
nest()%>%
mutate(test = map(data, ~prop.test(.$n, .$N)), # Z - test for propotion
sumr = map(test, broom::tidy))%>%
select(sumr)%>%
unnest()%>%
select(estimate1:p.value)%>%
rename(`H. schistosus` = estimate2,
`H. curtus` = estimate1)| H. curtus | H. schistosus | statistic | p.value |
|---|---|---|---|
| 0.4255319 | 0.179938 | 71.81006 | 0 |
H. curtus had as significatly higher mortality rate in bycatch than H. schistosus
Mortality rate by age class
Are juveniles more vulnerable to bycatch mortality than adults?
snakes%>%
filter(!is.na(Class))%>%
group_by(Species, Class)%>%
count(Condition.at.encounter..D.A.)%>%
#calculating mortality rate
mutate(N = sum(n),
prop.dead = n/N)%>%
filter(Condition.at.encounter..D.A. == "Dead")%>%
select(-Condition.at.encounter..D.A.)| Species | Class | n | N | prop.dead |
|---|---|---|---|---|
| Hydrophis curtus | Adult | 40 | 62 | 0.6451613 |
| Hydrophis curtus | Juvenile | 41 | 125 | 0.3280000 |
| Hydrophis schistosus | Adult | 138 | 648 | 0.2129630 |
| Hydrophis schistosus | Juvenile | 5 | 105 | 0.0476190 |
| Hydrophis schistosus | Neonate | 3 | 8 | 0.3750000 |
Testing difference in mortality rates across age classes within species
snakes%>%
filter(!is.na(Class))%>%
group_by(Species, Class)%>%
count(Condition.at.encounter..D.A.)%>%
# Caluclating mortality rate
mutate(N = sum(n),
prop.dead = n/N)%>%
filter(Condition.at.encounter..D.A. == "Dead")%>%
select(-Condition.at.encounter..D.A.)%>%
group_by(Species)%>%
nest()%>%
mutate(test = map(data, ~prop.test(.$n, .$N)), # Z - test for proportion
sumr = map(test, broom::tidy))%>%
select(sumr)%>%
unnest()%>%
select(estimate1:p.value)| Species | estimate1 | estimate2 | statistic | p.value |
|---|---|---|---|---|
| Hydrophis curtus | 0.6451613 | 0.328000 | 15.71186 | 0.0000738 |
| Hydrophis schistosus | 0.2129630 | 0.047619 | 17.68174 | 0.0001447 |
H. curtus had the highest mortality rate overall. H. schistosus juveniles had the lowest.
Mortality rate by sex
snakes%>%
filter(Sex != "")%>%
group_by(Species, Sex)%>%
count(Condition.at.encounter..D.A.)%>%
# Calculating mortality rate
mutate(N = sum(n),
prop.dead = n/N)%>%
filter(Condition.at.encounter..D.A. == "Dead")%>%
select(-Condition.at.encounter..D.A.)| Species | Sex | n | N | prop.dead |
|---|---|---|---|---|
| Hydrophis curtus | Female | 36 | 102 | 0.3529412 |
| Hydrophis curtus | Male | 52 | 110 | 0.4727273 |
| Hydrophis schistosus | Female | 62 | 304 | 0.2039474 |
| Hydrophis schistosus | Male | 72 | 314 | 0.2292994 |
Testing difference in mortality rates across sexes withing species
snakes%>%
filter(Sex != "")%>%
group_by(Species, Sex)%>%
count(Condition.at.encounter..D.A.)%>%
mutate(N = sum(n),
prop.dead = n/N)%>%
filter(Condition.at.encounter..D.A. == "Dead")%>%
select(-Condition.at.encounter..D.A.)%>%
group_by(Species)%>%
nest()%>%
mutate(test = map(data, ~prop.test(.$n, .$N)),
sumr = map(test, broom::tidy))%>%
select(sumr)%>%
unnest()%>%
select(estimate1:p.value)| Species | estimate1 | estimate2 | statistic | p.value |
|---|---|---|---|---|
| Hydrophis curtus | 0.3529412 | 0.4727273 | 2.6538718 | 0.1032980 |
| Hydrophis schistosus | 0.2039474 | 0.2292994 | 0.4448479 | 0.5047918 |
There was no difference in mortality rates across sexes in either species.
Mortality and female reproductive status
Are gravid females more susceptible to bycatch mortality than the rest of our sample?
HS%>%
filter(Sex == "Female",
Class == "Adult")%>%
group_by(Gravid)%>%
count(Condition.at.encounter..D.A.)%>%
# Calculating mortality rate
mutate(N = sum(n),
prop.dead = n/N)%>%
filter(Condition.at.encounter..D.A. == "Dead")%>%
select(-Condition.at.encounter..D.A.)| Gravid | n | N | prop.dead |
|---|---|---|---|
| 41 | 151 | 0.2715232 | |
| Yes | 16 | 85 | 0.1882353 |
Testing thing difference in mortality rate of H. schistosus females by reproductive state
HS%>%
filter(Sex == "Female",
Class == "Adult")%>%
group_by(Gravid)%>%
count(Condition.at.encounter..D.A.)%>%
mutate(N = sum(n),
prop.dead = n/N)%>%
filter(Condition.at.encounter..D.A. == "Dead")%>%
select(-Condition.at.encounter..D.A.)%>%
ungroup()%>%
nest()%>%
mutate(test = map(data, ~prop.test(.$n, .$N)),
sumr = map(test, broom::tidy))%>%
select(sumr)%>%
unnest()%>%
select(estimate1:p.value)| estimate1 | estimate2 | statistic | p.value |
|---|---|---|---|
| 0.2715232 | 0.1882353 | 1.629856 | 0.2017229 |
Observed reproductive cycle of H. schistosus
Proportion of gravid females in the sample
# percentage of gravid females in sample
HS%>%
count(Gravid)%>%
mutate(N = sum(n))%>%
mutate(prop.gravid = n/N)%>%
filter(Gravid == "Yes")%>%
select(-c(Gravid))| n | N | prop.gravid |
|---|---|---|
| 88 | 967 | 0.0910031 |
Number of gravid females encountered per year
# checking the number of gravid females per year
HS%>%
group_by(Year)%>%
filter(Gravid == "Yes")%>%
count(Gravid)%>%
spread(Gravid, n)%>%
rename(`n(Gravid Females)` = Yes)| Year | n(Gravid Females) |
|---|---|
| 2016 | 1 |
| 2018 | 75 |
| 2019 | 12 |
Monthly variation in proportion of gravid females
Proper sampling was only conducted in 2018/19 and hence only this period is used for analysis of reproductive cycles.
# calculating the proportion of gravid females per month
gravid <- HS%>%
mutate(Month = factor(Month, levels = month.name))%>%
filter(Year == "2018" | Year == "2019")%>% # only for 2018/19
group_by(Month)%>%
summarise(N = n(),
gravid = sum(Gravid == "Yes"),
prop.gravid = gravid/N)
gravid| Month | N | gravid | prop.gravid |
|---|---|---|---|
| January | 102 | 12 | 0.1176471 |
| February | 134 | 34 | 0.2537313 |
| March | 129 | 24 | 0.1860465 |
| April | 54 | 7 | 0.1296296 |
| May | 36 | 1 | 0.0277778 |
| October | 16 | 0 | 0.0000000 |
| November | 69 | 3 | 0.0434783 |
| December | 34 | 6 | 0.1764706 |
Describing change in the relative proportions of gravid females across months and years of sampling.
# plotting prop gravid per month
gravid%>%
mutate(Month = factor(Month, levels = month.name))%>%
complete(Month, fill = list(prop.gravid = 0))%>%
ggplot(aes(Month, prop.gravid))+
geom_point(size = 3)+
geom_path(aes(group = 1), size = 1, linetype = "dashed")+
geom_segment(aes(x = "April", xend = "April", y = 0 , yend = 0.02), #marking births
arrow = arrow(length = unit(0.25, "cm"), ends = "first"), size = 1)+
geom_text(aes(x = "April", y = 0.04, label = paste("Observed \n births")))+
geom_vline(aes(xintercept = "June"), size = 1)+#start of the monsoon ban
geom_vline(aes(xintercept = "August"), size = 1)+#end of the monsoon ban
geom_text(aes(x = "July", y = 0.20, label = "Monsoon ban"))+
scale_x_discrete(guide = guide_axis(n.dodge = 2))+
labs(y = "Proportion of gravid females")Pregnancy from Novemeber to May. Peak in Feb. Observed two live births in April.
Development of embryos and eggs
Sample size
embryos%>%
summarise(n.mothers = length(unique(Field.Code)),
n.embryos = length(unique(Embryo.Code)))| n.mothers | n.embryos |
|---|---|
| 29 | 235 |
Summary statistics of egg measurements
embryos%>%
select(Egg.Length..mm.., Egg.Width..mm.., Egg.Weigth..g.., Snout.to.Vent..cm., Embryo.Weight..g.)%>%
skimr::skim()%>%
skimr::yank("numeric")%>%
select(c(skim_variable, mean, sd))Variable type: numeric
| skim_variable | mean | sd |
|---|---|---|
| Egg.Length..mm.. | 39.71 | 11.72 |
| Egg.Width..mm.. | 32.23 | 9.01 |
| Egg.Weigth..g.. | 14.30 | 5.45 |
| Snout.to.Vent..cm. | 17.72 | 5.24 |
| Embryo.Weight..g. | 4.83 | 6.28 |
Sex ratios within clutches
embryos%>%
group_by(Field.Code)%>%
filter(Sex != "")%>%
summarise(prop.female = sum(Sex == "Female")/n())%>%
skimr::skim(prop.female)%>%
skimr::yank("numeric")%>%
select(c(skim_variable, mean, sd))Variable type: numeric
| skim_variable | mean | sd |
|---|---|---|
| prop.female | 0.53 | 0.31 |
Number of feamle embryos
embryos%>%
filter(Sex != "")%>%
summarise(females = sum(Sex == "Female"),
N = n())| females | N |
|---|---|
| 86 | 166 |
Testing clutch sex ratio
broom::tidy(prop.test(86, 166, p = 0.5))%>% # Z - test
select(estimate:conf.high)| estimate | statistic | p.value | parameter | conf.low | conf.high |
|---|---|---|---|---|---|
| 0.5180723 | 0.1506024 | 0.6979603 | 1 | 0.4395567 | 0.5957383 |
The sex ratio in clutches is not significantly different from 0.5.
Infertility rate
Percentage of eggs with no embryos present:
embryos%>%
count(Embryo)%>%
mutate(N = sum(n), rate = n*100/N)%>%
select(-N)| Embryo | n | rate |
|---|---|---|
| Absent | 5 | 2.12766 |
| Present | 230 | 97.87234 |
Matrotrophy
Do female H. schistosus expend energy in the development of embryos? or Are the eggs formed and yolk only provides nourishment?
Change in yolk weight with embryo weight
embryos%>%
ggplot(aes(Embryo.Weight..g., Yolk.Weight..g.))+
geom_point(size = 3)+
geom_smooth(method = "lm", linetype = "dashed", size = 1)+
labs(x = "Embryo weight (g)", y = "Yolk weight (g)")Linear model to test change in yolk weight with embryo development:
embryos%>%
select(Yolk.Weight..g., Embryo.Weight..g.)%>%
nest()%>%
# Yolk weight is normally distributed
mutate(mod = map(data, ~lm(Yolk.Weight..g. ~ Embryo.Weight..g., data = .)),
sumr = map(mod, broom::tidy),
stat = map(mod, broom::glance))%>%
select(sumr, stat)%>%
unnest()%>%
select(term:r.squared)| term | estimate | std.error | statistic | p.value | r.squared |
|---|---|---|---|---|---|
| (Intercept) | 13.0905565 | 0.2965229 | 44.14687 | 0 | 0.4071222 |
| Embryo.Weight..g. | -0.4526283 | 0.0450509 | -10.04705 | 0 | 0.4071222 |
Yolk weight decreases by 0.45 g for every 1g increase in embryo weight.
Change in egg mass with embryo development
embryos%>%
ggplot(aes(Embryo.Weight..g., Egg.Weigth..g..))+
geom_point(aes(col = Yolk.Weight..g.), size = 3)+
geom_smooth(method = "lm", linetype = "dashed", size = 1)+
scale_color_viridis(name = "Yolk Weight (g)")+
labs(x = "Embryo Weight (g)", y = "Total Egg Weight (g)")Linear model to test relation ship between total egg mass and embryo development:
embryos%>%
select(Egg.Weigth..g.., Embryo.Weight..g.)%>%
nest()%>%
# Total egg mass is normally distributed
mutate(mod = map(data, ~lm(Egg.Weigth..g.. ~ Embryo.Weight..g., data = .)),
sumr = map(mod, broom::tidy),
stat = map(mod, broom::glance))%>%
select(sumr, stat)%>%
unnest()%>%
select(term:r.squared)| term | estimate | std.error | statistic | p.value | r.squared |
|---|---|---|---|---|---|
| (Intercept) | 13.5441434 | 0.4114909 | 32.91481 | 0 | 0.2188866 |
| Embryo.Weight..g. | 0.3676922 | 0.0557912 | 6.59050 | 0 | 0.2188866 |
Total egg mass increases by 0.36 g for every 1g increase in embryo weight.
Note: Embryo weight is log-normally distributed.
Along with decrease of fat bodies of gravid females (Figure A1) as embryos develop indicates matrotrophic nutrition.
Reproductive effort
Does the amount of enery expended by female H. schistosus reduce with age?
# Cleaning reproductive effort data
re <- embryos%>%
left_join(snakes, by = c("Date", "Field.Code", "Species"))%>%
select(Date, Field.Code, Embryo.Code, Species, Egg.Length..mm..:Tail.Length..cm..x,
Egg.Weigth..g..:Sex.x, Snout.to.Vent..cm..y:Tail.Length..cm..y, Weight..g.,
-Body.Length..cm..x, - Body.Length..cm..y)%>%
filter(Species == "Hydrophis schistosus")%>%
rename(Embryo.SVL = Snout.to.Vent..cm..x,
Embryo.TL = Tail.Length..cm..x,
Embryo.Sex = Sex.x,
Female.SVL = Snout.to.Vent..cm..y,
Female.TL = Tail.Length..cm..y,
Mother.Wt = Weight..g.
)%>%
# Removing erroneous data point
filter(Female.SVL > 50)Increase in clutch size with female age
re_clutch <- re%>%
select(Field.Code, Egg.Weigth..g.., Mother.Wt, Female.SVL)%>%
group_by(Field.Code)%>%
summarise(clutch.size = n(),
Clutch.wt = sum(Egg.Weigth..g..),
Mother.wt = last(Mother.Wt),
Female.SVL = last(Female.SVL))%>%
mutate(rcm = Clutch.wt/(Mother.wt - Clutch.wt))
re_clutch%>%
ggplot(aes(Female.SVL, clutch.size))+
geom_point(size =3)+
geom_smooth(method = 'glm', method.args = list(family = 'poisson'))+
labs(x = "Female SVL (cm)", y = "Clutch Size")Generalised linear model to test change in clutch size with female SVL:
Clutch size is poisson distributed with mean 8.1052632 and variance 10.3216374
mcfadden <- function(x){
r2 <- 1 - (x$deviance/x$null.deviance)
}
re_clutch%>%
select(Female.SVL, clutch.size)%>%
nest()%>%
mutate(mod = map(data, ~glm(clutch.size ~ Female.SVL, data = ., family = "poisson")),
sumr = map(mod, broom::tidy),
stat = map(mod, broom::glance),
r.squared = map_dbl(stat, ~mcfadden(.)))%>%
select(sumr, stat,r.squared)%>%
unnest()%>%
select(c(term:p.value, r.squared))| term | estimate | std.error | statistic | p.value | r.squared |
|---|---|---|---|---|---|
| (Intercept) | -0.3764929 | 0.7174214 | -0.5247862 | 0.5997318 | 0.560094 |
| Female.SVL | 0.0256926 | 0.0073002 | 3.5194451 | 0.0004325 | 0.560094 |
The number of eggs borne by females increase by 1.0253151 for every 1 cm increase in female SVL.
Change in overall reproductive effort with age
re_clutch%>%
ggplot(aes(Female.SVL, rcm))+
geom_point(aes(col = clutch.size), size = 3)+
geom_smooth(method = "lm", linetype = "dashed")+
scale_y_continuous(name = "Relative clutch mass")+
scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
scale_color_viridis(name = "Clutch size")+
theme(legend.text = element_text(size = 10))Linear model to determine relationship between relative clutch mass and female SVL:
re_clutch%>%
select(Female.SVL, rcm)%>%
nest()%>%
mutate(mod = map(data, ~lm(rcm ~ Female.SVL, data = .)),
sumr = map(mod, broom::tidy),
stat = map(mod, broom::glance))%>%
select(sumr, stat)%>%
unnest()%>%
select(term:r.squared)| term | estimate | std.error | statistic | p.value | r.squared |
|---|---|---|---|---|---|
| (Intercept) | 0.6886163 | 0.4077131 | 1.6889727 | 0.1170168 | 0.063695 |
| Female.SVL | -0.0038362 | 0.0042458 | -0.9035133 | 0.3840327 | 0.063695 |
The overall reproductive effort does not change with female age.
Change in reproductive effort per embryo with female age
Does the effort per embryo change with female age?
# calculating reproductive effort per embruo
re_embryo <- re%>%
select(Field.Code, Embryo.Code, Egg.Weigth..g.., Mother.Wt, Female.SVL, Embryo.Sex, Embryo.SVL)%>%
group_by(Field.Code)%>%
mutate(clutch.wt = sum(Egg.Weigth..g..),
Female.wt = Mother.Wt - clutch.wt)%>%
group_by(Field.Code, Embryo.Code)%>%
summarise(inv = Egg.Weigth..g../Female.wt,# investment per embryo
Female.SVL = last(Female.SVL),
Embryo.Sex = Embryo.Sex,
Embryo.SVL = Embryo.SVL)%>%
ungroup()
# using residuals to control for embryo development
emsvlinv <- lm(inv ~ Embryo.SVL, data = re_embryo)
# Plotting
re_embryo%>%
modelr::add_residuals(emsvlinv)%>%
ggplot(aes(Female.SVL, resid))+
geom_point(aes(col = Embryo.SVL), size = 3)+
geom_smooth(method = "lm", linetype = "dashed", size = 1)+
scale_y_continuous(name = "Relative egg mass (residuals)")+
scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
scale_color_viridis(name = "Embryo SVL (cm)")Linear model to detemine relationship between female SVL and investment per embryo:
We have controlled for the effect of embryo development.
re_embryo%>%
select(Female.SVL, inv, Embryo.SVL)%>%
nest()%>%
mutate(mod = map(data, ~lm(inv ~ Female.SVL + Embryo.SVL, data = .)),
sumr = map(mod, broom::tidy),
stat = map(mod, broom::glance))%>%
select(sumr, stat)%>%
unnest()%>%
select(term:p.value, r.squared)| term | estimate | std.error | statistic | p.value | r.squared |
|---|---|---|---|---|---|
| (Intercept) | 0.0858802 | 0.0260524 | 3.296444 | 0.0017188 | 0.4433341 |
| Female.SVL | -0.0007894 | 0.0002588 | -3.050593 | 0.0035109 | 0.4433341 |
| Embryo.SVL | 0.0023248 | 0.0004171 | 5.573820 | 0.0000008 | 0.4433341 |
The relative egg mass (controlled for embryo development) reduces with female SVL.
Difference in reproductive effort for male and female embryos
Does female investment differ by sex of the embryo?
re_embryo%>%
filter(Embryo.Sex == "Male" | Embryo.Sex == "Female")%>%
modelr::add_residuals(emsvlinv)%>%
ggplot(aes(Female.SVL, resid, shape = Embryo.Sex))+
geom_point(aes(col = Embryo.SVL), size = 3)+
geom_smooth(method = "lm", linetype = "dashed", size = 1)+
scale_y_continuous(name = "Relative egg mass (residuals)")+
scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
scale_color_viridis(name = "Embryo SVL (cm)")+
scale_shape_discrete(name = "Embryo Sex")Linear model to test for the difference in female investment by embryo sex:
We controlled for the effect of embryo development.
re_embryo%>%
filter(Embryo.Sex == "Male" | Embryo.Sex == "Female")%>%
select(Embryo.Sex, inv, Embryo.SVL)%>%
nest()%>%
mutate(mod = map(data, ~lm(inv ~ Embryo.Sex + Embryo.SVL, data = .)),
sumr = map(mod, broom::tidy),
stat = map(mod, broom::glance))%>%
select(sumr, stat)%>%
unnest()%>%
select(term:p.value, r.squared)| term | estimate | std.error | statistic | p.value | r.squared |
|---|---|---|---|---|---|
| (Intercept) | 0.0165286 | 0.0090610 | 1.824151 | 0.0735636 | 0.3859103 |
| Embryo.SexMale | -0.0083114 | 0.0045803 | -1.814613 | 0.0750384 | 0.3859103 |
| Embryo.SVL | 0.0022764 | 0.0004455 | 5.109887 | 0.0000042 | 0.3859103 |
Female investment does not differ significantly with embryo sex.